Source code download:     [ .R ]

0.0.1 Introduction

systemPipeR provides flexible utilities for designing, building and running automated end-to-end analysis workflows for a wide range of research applications, including next-generation sequencing (NGS) experiments, such as RNA-Seq, ChIP-Seq, VAR-Seq and Ribo-Seq (H Backman and Girke 2016). Important features include a uniform workflow interface across different data analysis applications, automated report generation, and support for running both R and command-line software, such as NGS aligners or peak/variant callers, on local computers or compute clusters (Figure 1). The latter supports interactive job submissions and batch submissions to queuing systems of clusters.

systemPipeR has been designed to improve the reproducibility of large-scale data analysis projects while substantially reducing the time it takes to analyze complex omics data sets. Its unique features allow the creation of a uniform workflow interface and management system that allows the users to run selected steps, customize, and design entirely new workflows. Also, the package features take advantage of central community S4 classes of the Bioconductor ecosystem and command-line-based software support.

The main motivation and advantages of using systemPipeR for complex data analysis tasks are:

  1. Facilitates the design of complex workflows involving multiple R/Bioconductor packages
  2. Common workflow interface for different applications
  3. Makes analysis with Bioconductor utilities more accessible to new users
  4. Simplifies usage of command-line software from within R
  5. Reduces the complexity of using compute clusters for R and command-line software
  6. Accelerates runtime of workflows via parallelization on computer systems with multiple CPU cores and/or multiple compute nodes
  7. Improves reproducibility by automating analyses and generation of analysis reports

Figure 1: Relevant features in systemPipeR. Workflow design concepts are illustrated under (A). Examples of systemPipeR's visualization functionalities are given under (B).

A central concept for designing workflows within the systemPipeR environment is the use of workflow management containers. Workflow management containers allow the automation of design, build, run and scale different steps and tools in data analysis. systemPipeR adopted the widely used community standard Common Workflow Language (CWL) (Amstutz et al. 2016) for describing parameters analysis workflows in a generic and reproducible manner, introducing SYSargsList workflow control class (see Figure 2). Using this community standard in systemPipeR has many advantages. For instance, the integration of CWL allows running systemPipeR workflows from a single specification instance either entirely from within R, from various command-line wrappers (e.g., cwl-runner) or from other languages (, e.g., Bash or Python). systemPipeR includes support for both command-line and R/Bioconductor software as well as resources for containerization, parallel evaluations on computer clusters along with the automated generation of interactive analysis reports.

Figure 2: Overview of systemPipeR workflows management instances. A) A typical analysis workflow requires multiple single software (red), as well the description of the input data (green), and the expected outfiles and reports analysis (purple). B) systemPipeR provides multiple utilities to design and build a workflow, allowing multi-instance, integration of R code and command-line software, a simple and efficient annotation system, that allows automatic control of the input and output data, and multiple resources to manage the entire workflow. c) The execution of the analysis is independent of the design and build, enabling portability and more meaningful code sharing. systemPipeR provides options to execute a single step or multi-steps in a compute session environment, enabling scalability, resources to re-run one or multi-steps, checkpoints, and generation of execution reports that can track parameters and versions, providing transparency and data provenance.

An important feature of systemPipeR's CWL interface is that it provides two options to run command-line tools and workflows based on CWL. First, one can run CWL in its native way via an R-based wrapper utility for cwl-runner or cwl-tools (CWL-based approach). Second, one can run workflows using CWL’s command-line and workflow instructions from within R (R-based approach). In the latter case the same CWL workflow definition files (e.g. *.cwl and *.yml) are used but rendered and executed entirely with R functions defined by systemPipeR, and thus use CWL mainly as a command-line and workflow definition format rather than execution software to run workflows. In this regard systemPipeR also provides several convenience functions that are useful for designing and debugging workflows, such as a command-line rendering function to retrieve the exact command-line strings for each data set and processing step prior to running a command-line.

This overview introduces the design of a new CWL S4 class in systemPipeR, as well as the custom command-line interface, combined with the overview of all the common analysis steps of NGS experiments.

0.0.2 Workflow Management with SYSargsList

systemPipeR allows creation (multi-step analyses) and execution of workflow entirely for R, with control, flexibility, and scalability of all processes. Furthermore, the workflow execution can be integrated with compute clusters from R, accelerating results acquisition.

The flexibility of systemPipeR's new interface workflow management class is the driving factor behind the use of as many steps necessary for the analysis as well as the connection between command-line- or R-based software. The connectivity among all workflow steps is achieved by the SYSargsList workflow management class.

SYSargsList S4 class is a list-like container where each instance stores all the input/output paths and parameter components required for a particular data analysis step (see Figure 3).

The SYSargsList constructor function will generate the instances, using as data input initial targets files, as well as two-parameter files (for details, see below). When running preconfigured workflows, the only input the user needs to provide is the initial targets file containing the paths to the input files (e.g., FASTQ) along with unique sample labels. Subsequent targets instances are created automatically, based on the connectivity establish between the steps. The parameters required for running command-line software is provided by the parameter (*.cwl and *.yml)) files described below.

The class store one or multiple steps, allowing central control for running, checking status, and monitor complex workflows from start to finish. This design enhances the systemPipeR workflow framework with a generalized, flexible, and robust design.

Figure 3: Workflow steps with input/output file operations are controlled by SYSargsList objects. Each SYSargs2 instance is constructed from one targets and two param files or LineWise instance for R code-based. The only input provided by the user is the initial targets file. Subsequent targets instances are created automatically, from the previous output files. Any number of predefined or custom workflow steps are supported. One or many SYSargsList objects are organized in an SYSargsList container.

0.1 Getting Started

0.1.1 Installation

The R software for running systemPipeR can be downloaded from CRAN. The systemPipeR environment can be installed from the R console using the BiocManager::install command. The associated data package systemPipeRdata can be installed the same way. The latter is a helper package for generating systemPipeR workflow environments with a single command containing all parameter files and sample data required to quickly test and run workflows.

if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install("systemPipeR")
BiocManager::install("systemPipeRdata")

Please note that if you desire to use a third-party command line tool, the particular tool and dependencies need to be installed and executable. See details.

0.1.2 Loading package and documentation

library("systemPipeR")  # Loads the package
library(help = "systemPipeR")  # Lists package info
vignette("systemPipeR")  # Opens vignette

0.1.3 Load sample data and workflow templates

The mini sample FASTQ files used by this overview vignette as well as the associated workflow reporting vignettes can be loaded via the systemPipeRdata package as shown below. The chosen data set SRP010938 obtains 18 paired-end (PE) read sets from Arabidposis thaliana (Howard et al. 2013). To minimize processing time during testing, each FASTQ file has been subsetted to 90,000-100,000 randomly sampled PE reads that map to the first 100,000 nucleotides of each chromosome of the A. thalina genome. The corresponding reference genome sequence (FASTA) and its GFF annotation files (provided in the same download) have been truncated accordingly. This way the entire test sample data set requires less than 200MB disk storage space. A PE read set has been chosen for this test data set for flexibility, because it can be used for testing both types of analysis routines requiring either SE (single-end) reads or PE reads.

The following generates a fully populated systemPipeR workflow environment (here for RNA-Seq) in the current working directory of an R session. At this time the package includes workflow templates for RNA-Seq, ChIP-Seq, VAR-Seq, and Ribo-Seq. Templates for additional NGS applications will be provided in the future.

sytemPipeRdata::genWorkenvir(workflow = "rnaseq")
setwd("rnaseq")

0.1.4 Project structure

The working environment of the sample data loaded in the previous step contains the following pre-configured directory structure (Figure 4). Directory names are indicated in green. Users can change this structure as needed, but need to adjust the code in their workflows accordingly.

  • workflow/ (e.g. rnaseq/)
    • This is the root directory of the R session running the workflow.
    • Run script ( *.Rmd) and sample annotation (targets.txt) files are located here.
    • Note, this directory can have any name (e.g. rnaseq, varseq). Changing its name does not require any modifications in the run script(s).
    • Important subdirectories:
      • param/
        • param/cwl/: This subdirectory stores all the CWL parameter files. To organize workflows, each can have its own subdirectory, where all CWL param and input.yml files need to be in the same subdirectory.
      • data/
        • FASTQ files
        • FASTA file of reference (e.g. reference genome)
        • Annotation files
        • etc.
      • results/
        • Analysis results are usually written to this directory, including: alignment, variant and peak files (BAM, VCF, BED); tabular result files; and image/plot files.
        • Note, the user has the option to organize results files for a given sample and analysis step in a separate subdirectory.

Figure 4: systemPipeR’s preconfigured directory structure.

The following parameter files are included in each workflow template:

  1. targets.txt: initial one provided by user; downstream targets_*.txt files are generated automatically
  2. *.param/cwl: defines parameter for input/output file operations, e.g.:
    • hisat2/hisat2-mapping-se.cwl
    • hisat2/hisat2-mapping-se.yml
  3. Configuration files for computer cluster environments (skip on single machines):
    • .batchtools.conf.R: defines the type of scheduler for batchtools pointing to template file of cluster, and located in user’s home directory
    • *.tmpl: specifies parameters of scheduler used by a system, e.g. Torque, SGE, Slurm, etc.

0.1.5 Structure of targets file

The targets file defines all input files (e.g. FASTQ, BAM, BCF) and sample comparisons of an analysis workflow. The following shows the format of a sample targets file included in the package. It also can be viewed and downloaded from systemPipeR’s GitHub repository here. In a target file with a single type of input files, here FASTQ files of single-end (SE) reads, the first three columns are mandatory including their column names, while it is four mandatory columns for FASTQ files of PE reads. All subsequent columns are optional and any number of additional columns can be added as needed. The columns in targets files are expected to be tab separated (TSV format). The SampleName column contains usually short labels for referencing samples (here FASTQ files) across many workflow steps (e.g. plots and column titles). Importantly, the labels used in the SampleName column need to be unique, while technical or biological replicates are indicated by duplicated values under the Factor column. For readability and transparency, it is useful to use here a short, consistent and informative syntax for naming samples and replicates. To avoid problems with other packages or external software, it is recommended to use the basic naming rules for R objects and their components as outlined here. This is important since the values used under the SampleName and Factor columns are intended to be used as labels for naming columns or plotting features in downstream analysis steps.

Users should note here, the usage of targets files is optional when using systemPipeR’s new CWL interface. They can be replaced by a standard YAML input file used by CWL. Since for organizing experimental variables targets files are extremely useful and user-friendly. Thus, we encourage users to keep using them.

0.1.5.1 Structure of targets file for single-end (SE) samples

library(systemPipeR)
targetspath <- system.file("extdata", "targets.txt", package = "systemPipeR")
showDF(read.delim(targetspath, comment.char = "#"))

To work with custom data, users need to generate a targets file containing the paths to their own FASTQ files and then provide under targetspath the path to the corresponding targets file.

0.1.5.2 Structure of targets file for paired-end (PE) samples

For paired-end (PE) samples, the structure of the targets file is similar, where users need to provide two FASTQ path columns: FileName1 and FileName2 with the paths to the PE FASTQ files.

targetspath <- system.file("extdata", "targetsPE.txt", package = "systemPipeR")
showDF(read.delim(targetspath, comment.char = "#"))

0.1.5.3 Sample comparisons

Sample comparisons are defined in the header lines of the targets file starting with ‘# <CMP>.’

readLines(targetspath)[1:4]
## [1] "# Project ID: Arabidopsis - Pseudomonas alternative splicing study (SRA: SRP010938; PMID: 24098335)"                                                                              
## [2] "# The following line(s) allow to specify the contrasts needed for comparative analyses, such as DEG identification. All possible comparisons can be specified with 'CMPset: ALL'."
## [3] "# <CMP> CMPset1: M1-A1, M1-V1, A1-V1, M6-A6, M6-V6, A6-V6, M12-A12, M12-V12, A12-V12"                                                                                             
## [4] "# <CMP> CMPset2: ALL"

The function readComp imports the comparison information and stores it in a list. Alternatively, readComp can obtain the comparison information from the corresponding SYSargsList object (see below). Note, these header lines are optional. They are mainly useful for controlling comparative analyses according to certain biological expectations, such as identifying differentially expressed genes in RNA-Seq experiments based on simple pair-wise comparisons.

readComp(file = targetspath, format = "vector", delim = "-")
## $CMPset1
## [1] "M1-A1"   "M1-V1"   "A1-V1"   "M6-A6"   "M6-V6"   "A6-V6"   "M12-A12"
## [8] "M12-V12" "A12-V12"
## 
## $CMPset2
##  [1] "M1-A1"   "M1-V1"   "M1-M6"   "M1-A6"   "M1-V6"   "M1-M12"  "M1-A12" 
##  [8] "M1-V12"  "A1-V1"   "A1-M6"   "A1-A6"   "A1-V6"   "A1-M12"  "A1-A12" 
## [15] "A1-V12"  "V1-M6"   "V1-A6"   "V1-V6"   "V1-M12"  "V1-A12"  "V1-V12" 
## [22] "M6-A6"   "M6-V6"   "M6-M12"  "M6-A12"  "M6-V12"  "A6-V6"   "A6-M12" 
## [29] "A6-A12"  "A6-V12"  "V6-M12"  "V6-A12"  "V6-V12"  "M12-A12" "M12-V12"
## [36] "A12-V12"

0.2 Project initialization

To create a workflow within systemPipeR, we can start by defining an empty container and checking the directory structure:

sal <- SPRproject()
## Creating directory '/home/dcassol/z_tests/GEN242_toPush/content/en/tutorials/systempiper/rnaseq/.SPRproject'
## Creating file '/home/dcassol/z_tests/GEN242_toPush/content/en/tutorials/systempiper/rnaseq/.SPRproject/SYSargsList.yml'

Internally, SPRproject function will create a hidden folder called .SPRproject, by default, to store all the log files.

In this stage, the object sal is a empty container, except for the project information. The project information can be accessed by the projectInfo method:

sal
## Instance of 'SYSargsList': 
##  No workflow steps added
projectInfo(sal)
## $project
## [1] "/home/dcassol/z_tests/GEN242_toPush/content/en/tutorials/systempiper/rnaseq"
## 
## $data
## [1] "data"
## 
## $param
## [1] "param"
## 
## $results
## [1] "results"
## 
## $logsDir
## [1] ".SPRproject"
## 
## $sysargslist
## [1] ".SPRproject/SYSargsList.yml"

Also, the length function will return how many steps this workflow contains, and in this case, it is empty, as follow:

length(sal)
## [1] 0

0.3 Workflow Design

systemPipeR workflows can be designed and built from start to finish with a single command, importing from an R Markdown file or stepwise in interactive mode from the R console.

In the next section, we will demonstrate how to build the workflow in an interactive mode, and in the following section, we will show how to build from an R Markdown file.

New workflows are constructed, or existing ones modified, by connecting each step via appendStep method. Each SYSargsList instance contains instructions needed for processing a set of input files with a specific command-line and the paths to the corresponding outfiles generated.

The constructor function Linewise is used to build the R code-based step.

0.3.1 Build workflow interactive

This tutorial will demonstrate how to create a step in the workflow for running the short read aligner HISAT2 (Kim, Langmead, and Salzberg 2015).

The constructor function renders the proper command-line strings for each sample and software tool, appending a new step in the SYSargsList object defined in the previous step. For that, the SYSargsList function requires data from three input files:

- CWL command-line specification file (`wf_file` argument);
- Input variables (`input_file` argument);
- Targets file (`targets` argument).

In CWL, files with the extension .cwl define the parameters of a chosen command-line step or workflow, while files with the extension .yml define the input variables of command-line steps. Note, input variables provided by a targetsPE file can be passed on to a SYSargsList instance via the inputvars argument of the SYSargsList function.

appendStep(sal) <- SYSargsList(step_name = "hisat2_mapping", dir = TRUE, targets = "targetsPE.txt",
    wf_file = "workflow-hisat2/workflow_hisat2-pe.cwl", input_file = "workflow-hisat2/workflow_hisat2-pe.yml",
    dir_path = "param/cwl", inputvars = c(FileName1 = "_FASTQ_PATH1_", FileName2 = "_FASTQ_PATH2_",
        SampleName = "_SampleName_"))

For a overview of the workflow, we can check the object as follows:

sal
## Instance of 'SYSargsList': 
##     WF Steps:
##        1. hisat2_mapping --> Status: Pending 
##            Total Files: 72 | Existing: 0 | Missing: 72 
##          1.1. hisat2
##              cmdlist: 18 | Pending: 18
##          1.2. samtools-view
##              cmdlist: 18 | Pending: 18
##          1.3. samtools-sort
##              cmdlist: 18 | Pending: 18
##          1.4. samtools-index
##              cmdlist: 18 | Pending: 18
## 

Note that the workflow status is Pending, which means the workflow object is rendered in R; however, we did not execute the workflow yet.

Several accessor methods are available to explore the SYSargsList object.

Of particular interest is the cmdlist() method. It constructs the system commands for running command-line software as specified by a given .cwl file combined with the paths to the input samples (e.g. FASTQ files) provided by a targets file. The example below shows the cmdlist() output for running HISAT2 on the first PE read sample. Evaluating the output of cmdlist() can be very helpful for designing and debugging .cwl files of new command-line software or changing the parameter settings of existing ones.

For more details about the command-line rendered for each sample, it can be checked as follows:

cmdlist(sal, step = "hisat2_mapping", targets = 1)
## $hisat2_mapping
## $hisat2_mapping$M1A
## $hisat2_mapping$M1A$hisat2
## [1] "hisat2 -S ./results/M1A.sam  -x ./data/tair10.fasta  -k 1  --min-intronlen 30  --max-intronlen 3000  -1 ./data/SRR446027_1.fastq.gz -2 ./data/SRR446027_2.fastq.gz --threads 4"
## 
## $hisat2_mapping$M1A$`samtools-view`
## [1] "samtools view -bS -o ./results/M1A.bam  ./results/M1A.sam "
## 
## $hisat2_mapping$M1A$`samtools-sort`
## [1] "samtools sort -o ./results/M1A.sorted.bam  ./results/M1A.bam  -@ 4"
## 
## $hisat2_mapping$M1A$`samtools-index`
## [1] "samtools index -b results/M1A.sorted.bam  results/M1A.sorted.bam.bai  ./results/M1A.sorted.bam "

The outfiles components of SYSargsList define the expected output files for each step in the workflow; some of which are the input for the next workflow step, here next SYSargsList instance (see Figure 3).

outfiles(sal)
## $hisat2_mapping
## DataFrame with 18 rows and 4 columns
##              hisat2_sam       samtools_bam      samtools_sort_bam
##             <character>        <character>            <character>
## M1A   ./results/M1A.sam  ./results/M1A.bam ./results/M1A.sorted..
## M1B   ./results/M1B.sam  ./results/M1B.bam ./results/M1B.sorted..
## A1A   ./results/A1A.sam  ./results/A1A.bam ./results/A1A.sorted..
## A1B   ./results/A1B.sam  ./results/A1B.bam ./results/A1B.sorted..
## V1A   ./results/V1A.sam  ./results/V1A.bam ./results/V1A.sorted..
## ...                 ...                ...                    ...
## M12B ./results/M12B.sam ./results/M12B.bam ./results/M12B.sorte..
## A12A ./results/A12A.sam ./results/A12A.bam ./results/A12A.sorte..
## A12B ./results/A12B.sam ./results/A12B.bam ./results/A12B.sorte..
## V12A ./results/V12A.sam ./results/V12A.bam ./results/V12A.sorte..
## V12B ./results/V12B.sam ./results/V12B.bam ./results/V12B.sorte..
##              samtools_index
##                 <character>
## M1A  ./results/M1A.sorted..
## M1B  ./results/M1B.sorted..
## A1A  ./results/A1A.sorted..
## A1B  ./results/A1B.sorted..
## V1A  ./results/V1A.sorted..
## ...                     ...
## M12B ./results/M12B.sorte..
## A12A ./results/A12A.sorte..
## A12B ./results/A12B.sorte..
## V12A ./results/V12A.sorte..
## V12B ./results/V12B.sorte..

In an ‘R-centric’ rather than a ‘CWL-centric’ workflow design the connectivity among workflow steps is established by creating the downstream targets files automatically (see Figure 3). Each step can uses the previous step outfiles as an input, establishing connectivity among the steps in the workflow. By chaining several SYSargsList steps together one can construct complex workflows involving many sample-level input/output file operations with any combination of command-line or R-based software. Also, systemPipeR provides features to automatically and systematically build this connection, providing security that all the samples will be managed efficiently and reproducibly.

Alternatively, a CWL-centric workflow design can be used that defines all/most workflow steps with CWL workflow and parameter files. Due to time and space restrictions, the CWL-centric approach is not covered by this tutorial.

The following step generate data frame containing important read alignment statistics such as the total number of reads in the FASTQ files, the number of total alignments, as well as the number of primary alignments in the corresponding BAM files.

This constructor function LineWise requires the step_name and the R code-based under the code argument. The R code should be enclosed by braces ({}) and separated by a new line.

appendStep(sal) <- LineWise(code = {
    fqpaths <- getColumn(sal, step = "hisat2_mapping", "targetsWF", column = "FileName1")
    bampaths <- getColumn(sal, step = "hisat2_mapping", "outfiles", column = "samtools_sort_bam")
    read_statsDF <- alignStats(args = bampaths, fqpaths = fqpaths, pairEnd = TRUE)
    write.table(read_statsDF, "results/alignStats.xls", row.names = FALSE, quote = FALSE,
        sep = "\t")
}, step_name = "align_stats", dependency = "hisat2_mapping")

Also, for printing and double-check the R code in the step, the codeLine method can be used:

codeLine(sal, "align_stats")
## align_stats
##     fqpaths <- getColumn(sal, step = "hisat2_mapping", "targetsWF", column = "FileName1")
##     bampaths <- getColumn(sal, step = "hisat2_mapping", "outfiles", column = "samtools_sort_bam")
##     read_statsDF <- alignStats(args = bampaths, fqpaths = fqpaths, pairEnd = TRUE)
##     write.table(read_statsDF, "results/alignStats.xls", row.names = FALSE, quote = FALSE, sep = "\t")

One interesting feature showed here is the getColumn method that allows extracting the information for a workflow instance. Those files can be used in an R code, as demonstrated below

getColumn(sal, step = "hisat2_mapping", "outfiles", column = "samtools_sort_bam")
##                         M1A                         M1B 
##  "./results/M1A.sorted.bam"  "./results/M1B.sorted.bam" 
##                         A1A                         A1B 
##  "./results/A1A.sorted.bam"  "./results/A1B.sorted.bam" 
##                         V1A                         V1B 
##  "./results/V1A.sorted.bam"  "./results/V1B.sorted.bam" 
##                         M6A                         M6B 
##  "./results/M6A.sorted.bam"  "./results/M6B.sorted.bam" 
##                         A6A                         A6B 
##  "./results/A6A.sorted.bam"  "./results/A6B.sorted.bam" 
##                         V6A                         V6B 
##  "./results/V6A.sorted.bam"  "./results/V6B.sorted.bam" 
##                        M12A                        M12B 
## "./results/M12A.sorted.bam" "./results/M12B.sorted.bam" 
##                        A12A                        A12B 
## "./results/A12A.sorted.bam" "./results/A12B.sorted.bam" 
##                        V12A                        V12B 
## "./results/V12A.sorted.bam" "./results/V12B.sorted.bam"

0.3.2 Build workflow from R Markdown

The workflow can be created by importing the steps from an R Markdown file. As demonstrated above, it is required to initialize the project with SPRproject function.

importWF function will scan and import all the R chunk from the R Markdown file and build all the workflow instances. Then, each R chuck in the file will be converted in a workflow step.

sal <- SPRproject()
## Creating directory '/home/dcassol/z_tests/GEN242_toPush/content/en/tutorials/systempiper/rnaseq/.SPRproject'
## Creating file '/home/dcassol/z_tests/GEN242_toPush/content/en/tutorials/systempiper/rnaseq/.SPRproject/SYSargsList.yml'
sal <- importWF(sal, file_path = "systemPipeRNAseq.Rmd")
## Reading Rmd file
## 
##  ---- Actions ----
## Ignore none-R chunks at line: 25
## Checking chunk eval values
## Checking chunk SPR option
## Ignore non-SPR chunks: 34, 43, 61, 96, 120, 187, 204, 297, 656, 684, 702, 710, 721
## Parse chunk code
## Now importing step 'load_SPR' 
## Now importing step 'preprocessing' 
## Now importing step 'trimming' 
## Now importing step 'fastq_report' 
## Now importing step 'hisat2_index' 
## Now importing step 'hisat2_mapping' 
## Now importing step 'align_stats' 
## Now importing step 'bam_IGV' 
## Now importing step 'create_db' 
## Now importing step 'read_counting' 
## Now importing step 'sample_tree' 
## Now importing step 'run_edger' 
## Now importing step 'custom_annot' 
## Now importing step 'filter_degs' 
## Now importing step 'venn_diagram' 
## Now importing step 'get_go_annot' 
## Now importing step 'go_enrich' 
## Now importing step 'go_plot' 
## Now importing step 'heatmap' 
## Now importing step 'sessionInfo'

Let’s explore the workflow to check the steps:

stepsWF(sal)
dependency(sal)
codeLine(sal)
targetsWF(sal)

0.3.2.1 Third-party software tools

Current, systemPipeR provides the param file templates for third-party software tools. A list is provided in the following table.

Tool Name Description Step
bwa BWA is a software package for mapping low-divergent sequences against a large reference genome, such as the human genome.  Alignment
Bowtie2 Bowtie 2 is an ultrafast and memory-efficient tool for aligning sequencing reads to long reference sequences. Alignment
FASTX-Toolkit FASTX-Toolkit is a collection of command line tools for Short-Reads FASTA/FASTQ files preprocessing. Read Preprocessing
TransRate Transrate is software for de-novo transcriptome assembly quality analysis. Quality
Gsnap GSNAP is a genomic short-read nucleotide alignment program. Alignment
Samtools Samtools is a suite of programs for interacting with high-throughput sequencing data. Post-processing
Trimmomatic Trimmomatic is a flexible read trimming tool for Illumina NGS data. Read Preprocessing
Rsubread Rsubread is a Bioconductor software package that provides high-performance alignment and read counting functions for RNA-seq reads. Alignment
Picard Picard is a set of command line tools for manipulating high-throughput sequencing (HTS) data and formats such as SAM/BAM/CRAM and VCF. Manipulating HTS data
Busco BUSCO assesses genome assembly and annotation completeness with Benchmarking Universal Single-Copy Orthologs. Quality
Hisat2 HISAT2 is a fast and sensitive alignment program for mapping NGS reads (both DNA and RNA) to reference genomes. Alignment
Tophat2 TopHat is a fast splice junction mapper for RNA-Seq reads. Alignment
GATK Variant Discovery in High-Throughput Sequencing Data. Variant Discovery
STAR STAR is an ultrafast universal RNA-seq aligner. Alignment
Trim_galore Trim Galore is a wrapper around Cutadapt and FastQC to consistently apply adapter and quality trimming to FastQ files. Read Preprocessing
TransDecoder TransDecoder identifies candidate coding regions within transcript sequences. Find Coding Regions
Trinity Trinity assembles transcript sequences from Illumina RNA-Seq data. denovo Transcriptome Assembly
Trinotate Trinotate is a comprehensive annotation suite designed for automatic functional annotation of transcriptomes. Transcriptome Functional Annotation
MACS2 MACS2 identifies transcription factor binding sites in ChIP-seq data. Peak calling
Kallisto kallisto is a program for quantifying abundances of transcripts from RNA-Seq data. Read counting
BCFtools BCFtools is a program for variant calling and manipulating files in the Variant Call Format (VCF) and its binary counterpart BCF. Variant Discovery
Bismark Bismark is a program to map bisulfite treated sequencing reads to a genome of interest and perform methylation calls in a single step. Bisulfite mapping
Fastqc FastQC is a quality control tool for high throughput sequence data. Quality
Blast BLAST finds regions of similarity between biological sequences. Blast

Remember, if you desire to run any of these tools, make sure to have the respective software installed on your system and configure in the PATH. You can check as follows:

tryCL(command = "grep")

0.4 How to run a Workflow

0.4.1 Setup and requirements

To go through this tutorial, you need the following software installed:

  • R (version >=4.1.2)
  • systemPipeR package (version >=2.0.8)
  • Hisat2 (version >= 2.1.0)

If you desire to build your pipeline with any different software, make sure to have the respective software installed and available in your PATH. To make sure if the configuration is correct, on test it with:

tryCL(command = "hisat2")  ## 'All set up, proceed!'

0.4.2 Running the workflow

For running the workflow, runWF function will execute all the command lines store in the workflow container.

sal <- runWF(sal)

This essential function allows the user to choose one or multiple steps to be executed using the steps argument. However, it is necessary to follow the workflow dependency graph. If a selected step depends on a previous step(s) that was not executed, the execution will fail.

sal <- runWF(sal, steps = c(1, 3))

Also, it allows forcing the execution of the steps, even if the status of the step is 'Success' and all the expected outfiles exists. Another feature of the runWF function is ignoring all the warnings and errors and running the workflow by the arguments warning.stop and error.stop, respectively.

sal <- runWF(sal, force = TRUE, warning.stop = FALSE, error.stop = FALSE)

0.4.2.1 Parallelization on clusters

This section of the tutorial provides an introduction to the usage of the systemPipeR features on a cluster.

Alternatively, the computation can be greatly accelerated by processing many files in parallel using several compute nodes of a cluster, where a scheduling/queuing system is used for load balancing.

The resources list object provides the number of independent parallel cluster processes defined under the Njobs element in the list. The following example will run 18 processes in parallel using each 4 CPU cores. If the resources available on a cluster allow running all 18 processes at the same time, then the shown sample submission will utilize in a total of 72 CPU cores.

Note, runWF can be used with most queueing systems as it is based on utilities from the batchtools package, which supports the use of template files (*.tmpl) for defining the run parameters of different schedulers. To run the following code, one needs to have both a conffile (see .batchtools.conf.R samples here) and a template file (see *.tmpl samples here) for the queueing available on a system. The following example uses the sample conffile and template files for the Slurm scheduler provided by this package.

The resources can be appended when the step is generated, or it is possible to add these resources later, as the following example using the addResources function:

resources <- list(conffile=".batchtools.conf.R",
                  template="batchtools.slurm.tmpl", 
                  Njobs=18, 
                  walltime=120, ## minutes
                  ntasks=1,
                  ncpus=4, 
                  memory=1024, ## Mb
                  partition = "short"
                  )
sal <- addResources(sal, c("hisat2_mapping"), resources = resources)
sal <- runWF(sal)

Note: The example is submitting the jog to short partition. If you desire to use a different partition, please adjust accordingly.

0.4.3 Visualize workflow

systemPipeR workflows instances can be visualized with the plotWF function.

This function will make a plot of selected workflow instance and the following information is displayed on the plot:

- Workflow structure (dependency graphs between different steps); 
- Workflow step status, *e.g.* `Success`, `Error`, `Pending`, `Warnings`; 
- Sample status and statistics; 
- Workflow timing: running duration time. 

If no argument is provided, the basic plot will automatically detect width, height, layout, plot method, branches, etc.

plotWF(sal)

0.4.4 Technical report

systemPipeR compiles all the workflow execution logs in one central location, making it easier to check any standard output (stdout) or standard error (stderr) for any command-line tools used on the workflow or the R code stdout. Also, the workflow plot is appended at the beginning of the report, making it easier to click on the respective step.

sal <- renderLogs(sal)

0.4.5 Scientific report

systemPipeR compiles all the workflow execution logs in one central location, making it easier to check any standard output (stdout) or standard error (stderr) for any command-line tools used on the workflow or the R code stdout. Also, the workflow plot is appended at the beginning of the report, making it easier to click on the respective step.

sal <- renderLogs(sal)

0.5 Workflow initialization with templates

Workflow templates are provided via systemPipeRdata and GitHub. Instances of these workflows can be created with a single command.

0.5.1 RNA-Seq sample

Load the RNA-Seq sample workflow into your current working directory.

library(systemPipeRdata)
genWorkenvir(workflow = "rnaseq")
setwd("rnaseq")

0.5.1.1 Create the workflow

This template provides some common steps for a RNAseq workflow. One can add, remove, modify workflow steps by operating on the sal object.

sal <- SPRproject()
sal <- importWF(sal, file_path = "systemPipeRNAseq.Rmd", verbose = FALSE)

Workflow includes following steps:

  1. Read preprocessing
    • Quality filtering (trimming)
    • FASTQ quality report
  2. Alignments: HISAT2 (or any other RNA-Seq aligner)
  3. Alignment stats
  4. Read counting
  5. Sample-wise correlation analysis
  6. Analysis of differentially expressed genes (DEGs)
  7. GO term enrichment analysis
  8. Gene-wise clustering

0.5.1.2 Run workflow

sal <- runWF(sal)

Workflow visualization

plotWF(sal)

Report generation

sal <- renderReport(sal)
sal <- renderLogs(sal)

0.5.2 ChIP-Seq sample

Load the ChIP-Seq sample workflow into your current working directory.

library(systemPipeRdata)
genWorkenvir(workflow = "chipseq")
setwd("chipseq")

Workflow includes following steps:

  1. Read preprocessing
    • Quality filtering (trimming)
    • FASTQ quality report
  2. Alignments: Bowtie2 or rsubread
  3. Alignment stats
  4. Peak calling: MACS2
  5. Peak annotation with genomic context
  6. Differential binding analysis
  7. GO term enrichment analysis
  8. Motif analysis

0.5.2.1 Create the workflow

This template provides some common steps for a ChIPseq workflow. One can add, remove, modify workflow steps by operating on the sal object.

sal <- SPRproject()
sal <- importWF(sal, file_path = "systemPipeChIPseq.Rmd", verbose = FALSE)

0.5.2.2 Run workflow

sal <- runWF(sal)

Workflow visualization

plotWF(sal)

Report generation

sal <- renderReport(sal)
sal <- renderLogs(sal)

0.5.3 VAR-Seq sample

Load the VAR-Seq sample workflow into your current working directory.

library(systemPipeRdata)
genWorkenvir(workflow = "varseq")
setwd("varseq")

Workflow includes following steps:

  1. Read preprocessing
    • Quality filtering (trimming)
    • FASTQ quality report
  2. Alignments: gsnap, bwa
  3. Variant calling: VariantTools, GATK, BCFtools
  4. Variant filtering: VariantTools and VariantAnnotation
  5. Variant annotation: VariantAnnotation
  6. Combine results from many samples
  7. Summary statistics of samples

0.5.3.1 Create the workflow

This template provides some common steps for a VARseq workflow. One can add, remove, modify workflow steps by operating on the sal object.

sal <- SPRproject()
sal <- importWF(sal, file_path = "systemPipeVARseq.Rmd", verbose = FALSE)

0.5.3.2 Run workflow

sal <- runWF(sal)

Workflow visualization

plotWF(sal)

Report generation

sal <- renderReport(sal)
sal <- renderLogs(sal)

0.5.4 Ribo-Seq sample

Load the Ribo-Seq sample workflow into your current working directory.

library(systemPipeRdata)
genWorkenvir(workflow = "riboseq")
setwd("riboseq")

Workflow includes following steps:

  1. Read preprocessing
    • Adaptor trimming and quality filtering
    • FASTQ quality report
  2. Alignments: HISAT2 (or any other RNA-Seq aligner)
  3. Alignment stats
  4. Compute read distribution across genomic features
  5. Adding custom features to workflow (e.g. uORFs)
  6. Genomic read coverage along transcripts
  7. Read counting
  8. Sample-wise correlation analysis
  9. Analysis of differentially expressed genes (DEGs)
  10. GO term enrichment analysis
  11. Gene-wise clustering
  12. Differential ribosome binding (translational efficiency)

0.5.4.1 Create the workflow

This template provides some common steps for a RIBOseq workflow. One can add, remove, modify workflow steps by operating on the sal object.

sal <- SPRproject()
sal <- importWF(sal, file_path = "systemPipeRIBOseq.Rmd", verbose = FALSE)

0.5.4.2 Run workflow

sal <- runWF(sal)

Workflow visualization

plotWF(sal, rstudio = TRUE)

Report generation

sal <- renderReport(sal)
sal <- renderLogs(sal)

0.6 Version information

Note: the most recent version of this tutorial can be found here.

sessionInfo()
## R Under development (unstable) (2021-10-25 r81105)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.4 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /home/dcassol/src/R-devel/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] magrittr_2.0.3              batchtools_0.9.15          
##  [3] ape_5.6-2                   ggplot2_3.3.5              
##  [5] systemPipeR_2.1.27          ShortRead_1.53.1           
##  [7] GenomicAlignments_1.31.2    SummarizedExperiment_1.25.3
##  [9] Biobase_2.55.2              MatrixGenerics_1.7.0       
## [11] matrixStats_0.62.0          BiocParallel_1.29.21       
## [13] Rsamtools_2.11.0            Biostrings_2.63.3          
## [15] XVector_0.35.0              GenomicRanges_1.47.6       
## [17] GenomeInfoDb_1.31.7         IRanges_2.29.1             
## [19] S4Vectors_0.33.17           BiocGenerics_0.41.2        
## [21] BiocStyle_2.23.1           
## 
## loaded via a namespace (and not attached):
##  [1] nlme_3.1-157            bitops_1.0-7            webshot_0.5.3          
##  [4] httr_1.4.2              RColorBrewer_1.1-3      progress_1.2.2         
##  [7] tools_4.2.0             backports_1.4.1         bslib_0.3.1            
## [10] DT_0.22                 utf8_1.2.2              R6_2.5.1               
## [13] DBI_1.1.2               colorspace_2.0-3        withr_2.5.0            
## [16] prettyunits_1.1.1       tidyselect_1.1.2        compiler_4.2.0         
## [19] rvest_1.0.2             cli_3.2.0               formatR_1.12           
## [22] xml2_1.3.3              DelayedArray_0.21.2     bookdown_0.26          
## [25] sass_0.4.1              scales_1.2.0            checkmate_2.1.0        
## [28] rappdirs_0.3.3          systemPipeRdata_1.99.10 systemfonts_1.0.4      
## [31] stringr_1.4.0           digest_0.6.29           svglite_2.1.0          
## [34] rmarkdown_2.13          jpeg_0.1-9              pkgconfig_2.0.3        
## [37] htmltools_0.5.2         fastmap_1.1.0           htmlwidgets_1.5.4      
## [40] rlang_1.0.2             rstudioapi_0.13         jquerylib_0.1.4        
## [43] generics_0.1.2          hwriter_1.3.2.1         jsonlite_1.8.0         
## [46] crosstalk_1.2.0         dplyr_1.0.8             RCurl_1.98-1.6         
## [49] kableExtra_1.3.4        GenomeInfoDbData_1.2.8  Matrix_1.4-1           
## [52] Rcpp_1.0.8.3            munsell_0.5.0           fansi_1.0.3            
## [55] lifecycle_1.0.1         stringi_1.7.6           yaml_2.3.5             
## [58] zlibbioc_1.41.0         grid_4.2.0              parallel_4.2.0         
## [61] crayon_1.5.1            lattice_0.20-45         hms_1.1.1              
## [64] knitr_1.38              pillar_1.7.0            base64url_1.4          
## [67] codetools_0.2-18        glue_1.6.2              evaluate_0.15          
## [70] latticeExtra_0.6-29     data.table_1.14.2       remotes_2.4.2          
## [73] BiocManager_1.30.16     png_0.1-7               vctrs_0.4.1            
## [76] gtable_0.3.0            purrr_0.3.4             assertthat_0.2.1       
## [79] xfun_0.30               viridisLite_0.4.0       tibble_3.1.6           
## [82] ellipsis_0.3.2          brew_1.0-7

0.7 Funding

This project is funded by NSF award ABI-1661152.

References

Amstutz, Peter, Michael R Crusoe, Nebojša Tijanić, Brad Chapman, John Chilton, Michael Heuer, Andrey Kartashov, et al. 2016. “Common Workflow Language, V1.0,” July. https://doi.org/10.6084/m9.figshare.3115156.v2.
H Backman, Tyler W, and Thomas Girke. 2016. systemPipeR: NGS workflow and report generation environment.” BMC Bioinformatics 17 (1): 388. https://doi.org/10.1186/s12859-016-1241-0.
Howard, Brian E, Qiwen Hu, Ahmet Can Babaoglu, Manan Chandra, Monica Borghi, Xiaoping Tan, Luyan He, et al. 2013. “High-Throughput RNA Sequencing of Pseudomonas-Infected Arabidopsis Reveals Hidden Transcriptome Complexity and Novel Splice Variants.” PLoS One 8 (10): e74183. https://doi.org/10.1371/journal.pone.0074183.
Kim, Daehwan, Ben Langmead, and Steven L Salzberg. 2015. HISAT: A Fast Spliced Aligner with Low Memory Requirements.” Nat. Methods 12 (4): 357–60.